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Black hole binaries on non-eccentric orbits form an important subclass of gravitational wave 
sources, but it is a non-trivial issue to construct numerical initial data with minimal initial eccen- 
tricity for numerical simulations. We compute post-Newtonian orbital parameters for quasi-spherical 
orbits using the method of Buonanno, Chen and Damour (2006) and examine the resulting eccen- 
tricity in numerical simulations. Four different methods are studied resulting from the choice of 
Taylor-expanded or effective-one-body Hamiltonians, and from two choices for the energy flux. For 
equal-mass, non-spinning binaries the approach succeeds in obtaining low eccentricity numerical 
initial data with an eccentricity of about e = 0.002 for rather small initial separations of D > 10M. 
The eccentricity increases for unequal masses and for spinning black holes, but remains smaller than 
that obtained from previous post-Newtonian approaches. The effective-one-body Hamiltonian offers 
advantages for decreasing initial separation as expected, but in the context of this study also per- 
forms significantly better than the Taylor-expanded Hamiltonian for binaries with spin. For mass 
ratio 4:1 and vanishing spin, the eccentricity reaches e = 0.004. For mass ratio 1:1 and aligned 
spins of size 0.85M 2 the eccentricity is about e = 0.07 for the Taylor method and e = 0.014 for the 
effective-one-body method. 

PACS numbers: 04.25. D-, 04.25.dg, 04.25.Nx 



I. INTRODUCTION 

Gravitational wave astronomy holds the promise to 
open up an entirely new window into the universe, and 
a number of ground-based gravitational wave detectors 
are now collecting data in the quest for first direct detec- 
tion [IHI] . Binary black hole (BBH) mergers are expected 
to be primary sources of gravitational waves. Theoretical 
predictions for waves from BBH mergers are now becom- 
ing available for various merger scenarios based on nu- 
merical simulations as well as post-Newtonian (PN) ap- 
proximation methods, with increasing quality and range 
of validity (e.g. [6]). Among the likely astrophysi- 
cal scenarios is that of two black holes merging at the 
end-point of a non-eccentric, quasi-circular inspiral. Al- 
though eccentric, non-circular orbits are certainly possi- 
ble and of interest as well, non-eccentric orbits are ex- 
pected to be common since the emission of gravitational 
waves reduces eccentricity on time scales that are short 
compared to the life time of the binary. 

In order to predict waveforms for non-eccentric inspi- 
rals, numerical simulations for the Einstein equations 
have to face the issue of how to construct black hole ini- 
tial data with minimal initial eccentricity. This is a non- 
trivial task for several reasons. The main reason is that 
the initial data have to satisfy the constraint equations 
of general relativity. Typical initial data constructions 
allow us to specify the "bare" parameters of a multi- 
ple black hole system (masses, positions, momenta, and 
spins) , and the field content (the metric and extrinsic cur- 
vature on a hypersurface) is generated by solving elliptic 
equations. The solution process changes or dresses the 
bare parameters, e.g. the physical masses and the proper 



distance between the black holes are obtained as part of 
the solution. There are therefore two levels of indirect- 
ness in the initial data construction. Discrete parameters 
are translated into fields, and only after solving the con- 
straints the physical parameters are known. The goal is 
to obtain an accurate approximation to the physical sit- 
uation at some given instant of time, which as a matter 
of principle can only be an approximation and cannot 
be exact because the evolution prior to this time is not 
available. 

With regard to eccentricity, the question is how ini- 
tial parameters can be found that result in numerical 
binary inspirals with minimal eccentricity. In fact, most 
numerical simulations have to contend with a small but 
sometimes non-negligible amount of eccentricity, and it 
becomes a quantitative question to what extent the com- 
putation of faithful gravitational wave templates is ham- 
pered by unwanted eccentricity. For example, the Nu- 
merical Injection and Analysis (NINJA) project Tj could 
benefit from and could evaluate the importance of low- 
eccentricity initial configurations. 

Several methods are available to find parameters for 
approximately circular inspirals. A first approximation is 
to impose a type of quasi-circularity or quasi-equilibrium 
condition on the two black holes when solving the con- 
straints, which typically is implemented as an iteration 
of the orbital parameters and by repeatedly solving the 
constraints until the quasi-equilibrium condition is satis- 
fied within some given error (see e.g. [5]). However, the 
final judge of the eccentricity contained in initial data 
is to perform numerical evolutions and to measure the 
resulting eccentricity. This leads to the immediate sug- 
gestion to iterate initial parameters based on the eccen- 
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tricity observed in actual evolutions — which, of course, 
became an option only after numerical simulations of sev- 
eral orbits became possible. Ref. [5] suggests and applies 
a parameter iteration method based on evolutions with a 
pseudospectral code. The eccentricity is estimated after 
one or two orbits, and the run is restarted with improved 
parameters. The appeal of this method is its general- 
ity. It still has to be determined how well it can work 
with less accurate finite differencing codes. In particu- 
lar, due to artificial waves and gauge issues it is not clear 
how many orbits are required to measure the eccentricity 
with sufficient accuracy in the numerical simulation. For 
short simulations of only a few orbits, several iterations 
of the method may be computationally expensive. 

Another natural idea in this context is to obtain orbital 
parameters from PN methods. Once PN orbital param- 
eters for minimal eccentricity within the PN approxima- 
tion are found, they can be used as input for a constraint 
solution scheme within the full theory. Even simple im- 
plementations of this method can work surprisingly well. 
For example in [TU] (see also [TTJ [H]), it turns out that 
already a simple 3PN formula to compute the tangen- 
tial momentum P t for a PN circular orbit with vanishing 
radial momentum Pr, leads to approximately circular 
inspiral initial data when combined with the puncture 
method [12] to construct binary black hole initial data. 

Typical numerical simulations start shortly before the 
merger, say at an initial separation that allows 5 to 15 
orbits before merger. In this regime, the PN approach 
becomes increasingly inaccurate with decreasing separa- 
tion as the dynamic, strong field effects of the inspiral be- 
come severe. In order to improve on the PN parameters 
of [TU], we can look to more sophisticated PN approxi- 
mations, say involving the effective-one-body (EOB) ap- 
proach [14141 7j . In particular, the PN method should also 
supply a non- vanishing radial momentum Pr. Ironically, 
it is the use of "quasi-circular" initial parameters with 
vanishing radial momentum that leads to non-circular, 
eccentric inspirals, since by definition a spiral requires a 
non- vanishing radial momentum. 

We note in passing that there also exists a method to 
translate PN data including not only parameters P t and 
Pr, but also the PN field content to numerical initial 
data [TS]. The result is black hole puncture data with 
more realistic gravitational wave content [19] . a feature 
shared with the periodic standing wave approximation, 
e.g. [2D]. This feature is not essential for numerical evolu- 
tions, but in any case is not implemented by the methods 
discussed here. 

One way to obtain improved PN parameters for in- 
spirals with small eccentricity is to start with a simple 
initial guess and to perform evolutions of a few hundred 
orbits in the PN Hamiltonian formulation until the or- 
bit has circularized to a high degree due to the emission 
of gravitational waves [21] . In terms of computational 
effort, note that once the PN formulas for the Hamil- 
tonian system have been coded, the time integration is 
easily handled by a standard Mathematica function in 



less than a minute. The elegance of the PN evolution 
method is that for any reasonable initial guess the PN 
evolution leads to a circularized inspiralling orbit due to 
the physical process of wave emission, and this is achieved 
with little additional effort. 

For numerical simulations of equal mass binaries with- 
out spin, the residual eccentricity is reduced by up to 
a factor of five by the PN evolution method com- 
pared to simple 2PN quasi-circular parameters [22l [23] , 
with an eccentricity e < 0.002 at an initial separation of 
D = 11M. For comparison, several equal mass data sets 
are mentioned in [21], with e = 0.0016 the lowest value 
used for finite differencing codes, e = 0.008 a typical 
value prior to the recent improvements, but e < 5 x 10~ 5 
for the iteration method in a pseudospectral code [3]. 
For more general binaries with spin and unequal masses, 
there is some improvement of a typical initial guess using 
the PN evolution method, but the numerical eccentric- 
ity is not always as small as one may wish [55] [5S] . As 
suggested in [27] , it may be possible to add one iteration 
step to improve the situation for binaries with spin. 

The PN evolution method also demonstrates how large 
the eccentricity can be within the PN method itself if cer- 
tain simple prescriptions for quasi-circular orbit param- 
eters are used. On the other hand, the question arises 
how well a more sophisticated PN method can achieve 
orbital parameters for non-eccentric inspirals. In partic- 
ular, we want to ask the question whether there is a PN 
method that produces orbital parameters directly (with- 
out resorting to PN evolutions) that leads to reduced ec- 
centricity for numerical simulations. As explained above, 
we do not look for a method that reduces the eccentricity 
to zero, but any PN parameter proposal can be called suc- 
cessful e.g. if it reduces the eccentricity below the level of 
changes introduced by the constraint solving process. In 
any case, PN parameters can provide a good initial guess 
for an iteration based on full numerical simulations. 

One of the problems associated with the application of 
PN methods is that there is not one, but a multitude of 
approaches that differ in strategy as well as in small but 
important details. For example, two sets of PN initial 
parameters may be equivalent up to a certain PN order, 
but show differences in numerical simulations, which may 
be sensitive to the higher order terms that have been 
neglected. 

In this paper we focus on the proposal of Buonanno, 
Chen, and Damour [28 to compute quasi-spherical PN 
initial data. This work was based on the Hamiltonian for- 
mulation plus the known energy-loss rate and considered 
an adiabatic sequence of spherical orbits. That means 
the tangential momentum is determined by demanding 
a constant binary separation when using the conserva- 
tive part only. In a second step the radial component is 
derived by asking for the rate at which these spherical or- 
bits are (adiabatically) traversed by including the energy 
flux. The approach is capable of incorporating spin-orbit 
interactions to leading order, which may present a simple 
way of improving initial data for spinning binaries. 
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The present paper implements the suggestions and the 
algorithm of [28] and modifies it to be more directly appli- 
cable to dynamical variables instead of kinematical ones, 
i.e. ADMTT momenta instead of velocities. Following 
[2"B] we implement two variants of PN Hamiltonians, the 
usual Taylor-expanded and the EOB version. We also 
consider two versions of the energy flux, the classical 
Taylor-expanded flux for circularly orbiting bodies and 
a particular Pade approximant of it. When combined 
with the two choices of Hamiltonians this results in four 
different versions to study. 

In Section [II] we provide a detailed survey of the re- 
quired equations. We further examine the transition from 
EOB coordinates to ADMTT- type coordinates, which 
becomes necessary when extracting initial data from the 
EOB Hamiltonian. The initial data algorithm itself is 
described in Section IIIII Section IIVI is dedicated to the 
exploration of these data in numerical simulations with 
the BAM code. By measuring the orbital eccentricity 
obtained in numerical simulations we can argue which of 
the four data types is most appropriate. The investiga- 
tion includes mass ratios 1:1, 2:1, and 4:1 for vanishing 
spin, several cases of spin aligned to the orbital angular 
momentum for mass ratio 1:1, one anti-aligned case and 
one situation with arbitrarily chosen spins. Especially for 
binaries with spin the resulting eccentries are found to be 
quite large, but initial data involving the EOB Hamilto- 
nian still give reasonably small eccentricity. 



II. POST-NEWTONIAN EQUATIONS OF 
MOTION FOR SPINNING BINARY SYSTEMS 

This section provides the required equations of mo- 
tion of a black-hole binary system consisting of objects 
of masses m a , positions X a , momenta P a and spins S a 
(a = 1,2). For our purposes it is sufficient to restrict 
considerations to the center-of-mass dynamics, where 
P = P-y = — P 2 - Following [55] we use a Hamiltonian 
description in either ADMTT gauge or in EOB coordi- 
nates. 

The Hamiltonian formalism turns out to be very use- 
ful when working with the canonically conjugate position 
and momentum variables and in distinguishing conser- 
vative and radiative effects. The system's equations of 
motion take the form 



dx -ix m- dH 



dp 

~dt 
dS a 



= {P,H}+F = - 



dH 
dX ' 

dt = {*.,*} = ^xs B , 



F, 



(1) 
(2) 
(3) 



where F labels the non-conservative force and x denotes 
the usual vector cross product. The conservative part 
consists of an orbital and a spin contribution and we will 
use the total Hamiltonian formed by the sum of three 
components, 



H(X, P, St, S 2 ) =H°(X, P) + H so (X, P, Si, S 2 ) 

+ H SS (X,P,S 1 ,S 2 ). (4) 



The individual terms are presented in the subsequent 
paragraphs. 



Orbital Hamiltonian 



We consider two versions of the orbital contribu- 
tion H°, the standard 3PN accurate Taylor-expanded 
Hamiltonian (TH) derived by Damour, Jaranowski, and 
Schafer [151 l2"9H31| . and the effective-one-body Hamil- 
tonian (EH) given by Buonanno, Damour, Jaranowski, 
and Schafer at 3PN order in [T3 HH] (given first at 
2PN in [H]). It is crucial to note that the two Hamil- 
tonians are defined for different coordinate systems. 
The Taylor-expanded Hamiltonian uses ADMTT-type 
coordinates denoted by X, _P, where, however, a small 
deviation from the actual ADMTT gauge arises at the 
3PN- level [25] . The EOB coordinates, indicated by 
primes X',P', are related to the ADMTT coordinates 
by a canonical transformation. The character of this 
transformation will be investigated in Section ed D| below. 
Some equations are valid for both approaches, TH and 
EH, and in these cases the primes are omitted as in 
Eqns. 



Both versions of the Hamiltonian will be given ex- 
plicitly. By introducing reduced coordinates q = 
X/(GM) = (Xi - X 2 )/{GM) and p = P/fj,, where 
M = mi + fn 2 and fi = m\m 2 /M, the Taylor-expanded 
version reads 
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^Taylor(^P)=A^ 2 + M 



H N (q,p) + \H 1PN (q,p) + 4#2Pn(<?,p) + 4#3Pn(<?,p) 



with 



#n(<?,p) = 



p2 _ 1 

Y~q> 



i?ip N (g,p)=^(3^-i)(p 2 ) 2 -^ 



(3 + v)p 2 + v{npf 



1 1 

9 + 2?' 



(5) 

(6) 
(7) 



H 2PN (Q,P) =y & (l - 5i/ + 5z/ 2 ) (p 2 ) 3 + 1 



5 - 20i/ - 3^ 2 ) (p 2 ) 2 - 2v 2 {npfp 2 - 'dv 2 (np) 



H 3 PN(q,p) 



+ 



l 



3^(np) 2 + (5 + 8v)p 2 



q z 4 



^ (-5 + 35^ - 7(V 2 + 35^) (p 2 ) 4 + i 



(-7 + 42^ - 53z/ - 5z/ 3 )(p 2 ) 3 



+ (2 - 3^)z/ 2 (np) 2 (p 2 ) 2 + 3(1 - i/)i/ 2 (np)V - ^{npf 



-^(-27 + 136^ + 109i/ 2 )(p 2 ) 2 + -^(17 + 30i/)^(np) 2 p 2 + ^(5 + 43z/)i/(np) 4 
lb lb 12 



25 / 1 2 335 \ 23 , 
b4 48 



• ; I bi" i")""""}^ 



1 /109 21 o , 

- + 7T 2 V 

8 V 12 32 ! 



1 



(8) 



(9) 



where additionally q = \q\, n = q/q and v = n/M have been used. We already plugged in the appropriate values of 
the ambiguity parameters w sta tic an d wiunctic- 

The EOB Hamiltonian is given by the expression 



H^ OB {q',p') = Mc 2 x l + 2u 



H cf i(q',p') - /ic s 
fie 2 



(10) 



with 



H eS (q>,p>) =»c?JA(q>) 



c 2 \D(q') 

where A(q') is the Pade-resummed function 
A(q') 



' + ^ + ~ \) + ^ {Z1 iP ' 2)2 + Z2 P ' 2{n ' ■ P ' )2 + Z3 {n ' ■ P ' )4) 



q' 3 (8-2v) + jjq' 2 (a 4 + 8is - 16) 



g' 3 (8 - 2v) + ^q l2 {a A + Av) + ±q>{2a 4 + 8u) + ±4(a 4 + z/ 2 ) 

I 



(11) 



(12) 



and the remaining quantities are 
" '94 41 



As before we make use of reduced variables in the formu- 
lation of the EOB Hamiltonian. 



a 4 



7T 

3 32 



Zl 
V 



(13) 



^') = l-4 6 S + ^^ + ^ + (3,-2b) 
c 4 q' z c° L v v 



v 

7/3 : 



(14) 

zi = z 2 = and z 3 = 2(4 - Zv)v. (15) 
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B. Spin Hamiltonian 

The spin part will be considered to leading order only. 
These terms can be separated into spin-orbit and spin- 
spin interactions and have been known for a long time, 
see for example [28, 32J. We have 



Hs!S 2 
Hs 2 s 2 



= 1 



, G SeS ■ L 

3 m 2 

4 mi 



Si 



1 



3 mi 

4 m 2 



S 2 , 



(16) 
(17) 



"Hi Si 

GJ_ 
G 1 

1 ^ 2. 



+ H Sl s 2 +H S2 s 2 , (18) 
3(S 1 -N)(S 2 -N)-(S 1 -S 2 )], (19) 

3(Si-JV) (Si • JV) - (Si • Si) 



m 2 



mi 
(20) 

(21) 



Here, L = X x P is the orbital angular momentum 
and N is the same unit vector as n above. We assume 
the Newton- Wigner spin supplementary condition, which 
provides the notion of the spin vector and affects the def- 
inition of the bodies' worldlines. 

In what follows Eqns. ( 16 )-( 20 ) are presumed to hold 



for ADMTT as well as EOB coordinates. This, however, 
is an approximation as already noted in [28] . There exists 
a more sophisticated and potentially more accurate way 
of including spin effects in the EOB framework (T7] , but 
the Hamiltonian associated with this approach contains 
spin-orbit and spin-spin effects mixed in a complicated 
way, and it seems that these terms cannot be separated 
in a consistent manner. The initial-data algorithm we 
present, however, relics on neglecting the spin-spin terms 
for the sake of obtaining spherical orbits. Hence, in the 
EOB case we follow [28] by just employing the same sim- 
ple, additive spin Hamiltonian as in the ADM framework. 
We discuss the consequences of this approximation below 
when analyzing the canonical transformation between the 
two coordinate systems. 



C. Radiation-reaction force 

The radiation-reaction-force term F to be applied in 
Eq. ([2| is the expression derived by [58] . It has the fea- 
ture that it is formulated in terms of dynamical quanti- 
ties (X, P) instead of kinematical ones (X, X). It also 
includes spin interactions to linear order and reads 



1 dE 
~oj\L\ ~dt J 



15 U L 2 R 



61 + 48 rn l)p.S x 
mi J 



61 + 48 



mi 
m 2 



P-Sa >L. 



(22) 



The invariant velocity parameter 



(23) 



based on the orbital frequency u) = (p has been intro- 
duced, although this brings the kinematical quantities 
back into play. 

Eq. (22 1 is valid for (quasi-)circular orbits only, which 
is a result of simplifications in its derivation and due 
to the fact that the quantity dE/dt is the well-known 
energy- flux function of circularly orbiting masses [281 1331 - 
[37]. In particular, its standard form (up to 3.5PN order) 
does not account for any eccentricity parameter, 



dE 
dF 



+ \S4y) + ftssK + h(y)v* + h{v)v% (24) 



These equations are also available for non-vanishing ec- 
centricity in [35]. The expansion coefficients read 



h{r> 



1247 35 



336 12 



v. 



Air, 



44 711 9271 



65 



v H v 

9072 504 18 



-( 



8191 583 



V 672 24 



6 643 739 519 
69 854400 



7776 



16 2 


1712 




105 7E 




94403 






"IT J 


3024 



v 2 - 



775 
324' 



1712 
T05~' 



16 285 214 745 



504 



1728 



-v + 



193 385 
3024 



V 2 | 7T, 



fw 
/3SO 

/4SS 



(25) 

where 7e = 0.577215 ... is Euler's gamma and L gives 
the unit vector in the direction of L. Note that the term 



/ll 5to 2 \£-Si 




5 mi \ L ■ S 2 


V 4 + Ami) M 2 




4m 2 ) M 2 


V 


289{L-Si)(L 


s 2 )- 


103(Si • S 2 ) 


48 m\ m 2 


+ 0(S 2 i)- 
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/4SS is not fully known yet, but we will not be using it 
anyway. 



The force expression (22) actually does not enter the 



initial data algorithm presented in Section |III| below, 
however the energy flux (24) does. The fact that the 



flux is unable to describe eccentricity damping properly 
is not supposed to affect the algorithm. This is because 
we will merely ask for the shrinkage rate of circular (resp. 
spherical) orbits at one particular radius. In contrast, 
the application of the flux and the force F in PN sim- 
ulations with noticeable (e.g. spin induced) eccentricity 
might give inaccurate results. In such cases we may want 
to consider time dependent flux expressions or time av- 
erages over elliptical orbits [33] , 

Concerning the PN convergence of the flux function, 
several authors have shown that the Pade-resummation 
technique is capable of improving convergence, in par- 
ticular in the test-mass limit. A recent paper, however, 
claims that no advantages can be found using the Pade- 
resummation approach for several specific examples [35] . 
In order to examine this issue in the context of our eccen- 
tricity study, we implement the standard flux ( 24 ) , but 



also consider a direct P-approximant of the flux. The 
final expression is too lengthy to be shown here. Since 
computer-algebra software provides an efficient way to 
perform the resummation, we restrict ourselves to a brief 
summary of the required procedure. A more detailed de- 
scription can be found in [5UJ EEQ • 

Padc's method defines the resummation of a plain Tay- 
lor series. Referring to (24), we first separate the spin- 



dependent expressions and only resum the non-spinning 
part. In the remainder we find a logarithmic term which 
is factored out in a normalized form before the actual 
resummation (Eq. (38) in |41|). The normalization pa- 
rameter xlso can be taken from [10] . Eq. (3.23). Another 
factorization (Eqns. (39) and (44) in [41 ) takes care of a 
putative pole in the expansion near the light ring, whose 
position ngbt ring is given by Eq. (3.22) in [3D] to 2PN 
accuracy. Thus, finally, one factor of the flux function is 
a plain Taylor expansion and is resummed with a near- 
diagonal (lower-diagonal) Pade approximation. 

As a result we have a choice between two versions of 



dE/dt that are to be inserted in Eq. (22). We will refer 
to the different fluxes as Taylor flux (TF) and Pade flux 
(PF). 



D. Mapping EOB to ADMTT coordinates 

The link between ADM and EOB coordinates is given 
by a canonical transformation with generating function 

g(q,q r ), 



Pl dq* = p\ dq' 1 + dg{q, q') 



(26) 



where gid{q,q') = would correspond to the identity 
transformation [21 (TB] . Via a Legendre transformation 
we arrive at a more convenient generator G(q,p') = 



g(q,q') +p'i q H yielding 

ft d 9 l + g' l d^ = dG( 9 ,p')- (27) 
Decomposition into the identity map plus PN corrections, 
G{q, P ')=q l p' l + G{q,p'), (28) 

G{q,p r ) = -oGipnC^p') + 4g 2PN (<?,p') 

c z c 4 

+ ~G 3PN (q,p'), (29) 

was the key to the determination of the generating func- 
tion, which was accomplished by [T^JUFJ. Their result, 
which is not unique when including the 3PN level, reads 



Gipn(<7,p') ={q-p') aiP 



g 2 pn(<?,p') ={q-p') 



g 3 pn(q,p') ={q-p') 



II . a 2 



(30) 



hp' 4 + - (hp' 2 + b 3 (n ■ p') 



r\2\ 



(31) 



cip 



/(> 



-I- Ci(n- p' 



with 



Oi 

h 



+ ^ 

q 3 



v 
~2' 

-(zz + 3^), 

;(8^ + 3l/ 2 ) 

f 



1 (c2p' 4 + C3P' 2 (n-p') 2 
Q 

+ — (c 5 p /2 +c 6 (n-p') 2 ) 
q l 

(32) 



a 2 = l + -, 



(33) 



fc 2 = -(2^-5^ 2 ), 

W = 7(l-7i/ + ^ 2 ), (34) 



Cl 
C2 
C3 
C 4 
C5 
C6 
C7 



1G 
I 

16 
1 

24 
1 

16 
1 

16 
1 

"48 



(1 + Sis + 5is z )is, 
(l + 2i/- llv 2 )v, 
(12 + 48^ + 23^ 2 )^, 
(24 + 7is)is 2 , 
(13- Wis + 6is 2 )is, 
(115 + 116^ -26^ 2 )^, 



(35) 



1 2 155 

■7T + — I IS ■ 



64 



24 



1 



-is 



-is 



The sought-after transformation EOB — > ADM is 
<_ ,< dG{q,p') 



Pi=Pi + — Q^l — , (36) 
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which gives an implicit relation that can be solved nu- 
merically by a root-finding algorithm. Another possibil- 
ity to resolve the implicity is to repeatedly plug in the 
first equation in order to eliminate the unknown value q 
to sufficient PN order, yielding an explicit formula. By 
rearranging ( 36 ) we can calculate the inverse transforma- 



tion, ADM — > EOB, in an analogous manner. The inverse 
transformation is convenient when discussing properties 
of the transformation and the associated errors, for ex- 
ample we can start with ADM data from the class of 
circularly orbiting BH binaries and study the result on 
the EOB side. 

In order to illustrate the effect of the transformation, 
we consider the equal-mass case and choose, as a result 
from our initial data algorithm, the values 



q =(10,0,0), 

p = (-0.00394009,0.384516,0), 



(37) 



for a separation of D = 10M (with M = c = G = 1). 
This translates to 



q' = (10.9445,0.00318625,0), 
p' = (-0.00397814,0.351514,0). 



(38) 



when using the explicit version of ADM — > EOB at 3PN 
order. Note that the transformation mixes corresponding 
components of (q,p) and (q',p'), so the non- vanishing 
y-component of p causes a non-zero y-component of q' . 
The intermixture, however, is not a contradiction to the 
role of q and q' as position variables and p and p' as the 
conjugate momenta. 

For equal masses, Fig. [T] illustrates the actual relation 
between both coordinate types in the range from 6M to 
1AM, imposing non-eccentric momenta throughout. In 
particular, Fig. [l]a reveals that \q\ and \q'\ (resp. R and 
R') lie apart from each other by almost 1M. 

However, the validity of the canonical transforma- 
tion depends on the assumptions \q' — q\/\q\ <C 1 and 
\p' ~ pI/IpI ^ 1- The fulfillment of these conditions can 
be considered critical for the close binary separations un- 
der consideration, which might result in severe problems 
when using the spin and flux equations with EOB vari- 
ables. Namely, as we have just found there is a significant 
difference between the radial distances R and R' , as well 
as some deviation in the direction vectors TV, TV'. These 
quantities appear in Eqns. (fl6l-(20|. However, only (16) 



is incorporated in the initial data algorithm. Moreover, 
the transformation dp A dq = dp' A dq' corresponds to 
dH A dt = dH' A dt' and implies a change in the time 
scale [13] . This also affects the orbital frequency ui which 
is used within the flux function (24). The quantity L, 
however, is the same in both systems. The notion of 
spin vectors is supposed to be unaltered up to 1PN level, 
but modifications occur at higher PN orders. For a more 
detailed investigation we refer to the recently published 
reference |42j . where the spin vectors are made subject 
of an additional canonical transformation. 



Although we are aware of these issues, we note that ap- 
plying suitable corrections in order to avoid these incon- 
sistencies would be rather involved and would probably 
be unnecessarily sophisticated for our purpose. As al- 
ready discussed, further issues arise anyway when trans- 
lating the initial parameters obtained this way into the 
actual initial data (i.e. metric and curvature fields). In 
the end, the measured eccentricity decides whether our 
data - obtained under the specified simplifications - are 
appropriate. 



III. ALGORITHM FOR QUASI-SPHERICAL 
INITIAL DATA 



If we take into account only the conservative part of 
the orbital post-Newtonian dynamics, H°, together with 
the leading-order spin-orbit interaction term, H^o, it is 
possible to find solutions which feature a constant coor- 
dinate separation R. These solutions are called spher- 
ical orbits. In general, this would no longer be true if 
we additionally included spin-spin coupling via -f/sS) al- 
though there might exist very special spin constellations 
for which a spherical orbit can be found. In the following 
we only consider spin-orbit coupling (with one exception 
mentionend below). 

Imposing spin-orbit coupling we find two additional 
conserved quantities, the magnitude of the angular mo- 
mentum L = \L\ and 



XL 



GAP 



(39) 



which gives the projection of the spins onto the direction 
of L. Changing to polar coordinates with L pointing 
along the z-axis, the total Hamiltonian depends on only 
four quantities, 



n , G Gfvl L V r, 

H(R, P R , L, X l) = H°(R, P R , L + 2^ £± 

cr c R a 



where we used the relation 



P 2 = Pi + Pi = Pr+^ 



(40) 
(41) 



and P v = L. Note that only the z-components of the 
spins enter. Beyond these considerations for the conser- 
vative part, radiation reaction will be included separately. 

The challenge for the initial data algorithm is to find 
tangential and radial momentum of the objects for a 
given R and xl such that no radial oscillations occur 
during evolution with the dissipative PN equations of 
motion. 

Note that demanding a certain initial separation R is 
a gauge dependent statement, which, however, is con- 
venient for our purposes since the ADMTT gauge re- 
sembles the 1+log/gamma-freezing gauge used in stan- 
dard puncture evolutions. Under these circumstances, 
the algorithm takes a somewhat simpler form when us- 
ing R in order to specify the binary separation than 
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FIG. 1: (Color online) Character of the canonical transformation ADM — > EOB for spherical initial data. The dotted line 
indicates the identity map q' = q and p' — p. The plots on the right compare the explicit and implicit solution method of 
the transformation (36 1 and apparently there is convergence between the two versions with increasing PN order. However, 



convergence is not achieved conclusively, and the observed differences arc used to define an error bar for the transformation. 
Since the explicit transformation is an approximation (truncated at finite order) to the implicit relation, it is less accurate. 



when specifying the invariant orbital frequency oj in- 
stead, as originally done by Buonanno et al. [28 . Fur- 
thermore, it eases the specification of spins, which can 
now be set with respect to the dynamical quantity L = 
Xq x P /\Xq x Pol instead of the kinematical variable 
[£ N ] =I x X /\X x X |. We use the subscript "0" 
to stress that the initial value is meant. 

The algorithm itself consists of two steps: (i) Deter- 
mine [Pt]o as required for spherical orbits by considering 
the conservative motion only, (ii) Calculate the compo- 
nent [Pr]o by means of the energy flux in order to get a 
quasi-spherical inspiral. 



In the first step of the procedure we determine [Pt]o, 
or equivalently L — R ■ [i"t]o, by demanding constant 
radius, 



= [R] 



dH(R 7 P R ,L, XL ) 
dP R 



(42) 



Since Pr enters quadratically into H this is equivalent to 
[Pr]o - 0. (43) 

The value of \l is known and thus there is only one 
remaining degree of freedom, namely the value of Lq . It 
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is fixed by demanding 



= [Pr]q 



dH(R,P R ,L, X L) 
OR 



(44) 



a condition which has to be solved numerically. A start- 
ing value for an iteration procedure can be provided by 
the 3PN-accurate tangential momentum for circular or- 
bits of non-spinning objects, 



3 3PN 



1/2 



1 
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c 1 16 
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,n \ (GM\ 

4556)^ + 104i/ 2 



7/2 



(45) 



a formula given in [10] . 

The second step incorporates radiation reaction and 
consequently a non-vanishing R. Assuming adiabatic in- 
spiral, that is a sequence of spherical orbits whose shrink- 
age in radius is determined by the energy flux, one hnds 



[dE/dt] 
[{dE/dR) 



sphjo 



(46) 



where dE/dt is given by Eq. (24) or its Pade resumma- 
tion. Here the knowledge of the initial orbital frequency 
uio is required and can be provided by the evolution equa- 
tion 



w 



'dH ' 




dH' 










(47) 



The initial conditions used in the RHS of this equation, 
i.e. the components of the momentum, are taken from 
the first step, in particular [Pr]o = 0. The other quan- 
tity appearing in ( 46 1 is the energy difference between 



neighboring spherical orbits and can be calculated as 

dH fdP R \ 



(-) 



dH 
~dR 



dR 

d~R, 



jph, 



dH 
~d~L 



dL 
dR 



dPR V dR J sph 

sph 9 xl v di? ; sph 



The first two terms vanish, because for spherical or- 
bits dH/dR — —Pr — holds and consequently 
dPji = 0, too. The other two summands are problem- 
atic in that there appear the unknowns (dL/di?) sp h and 
(dxi/di?) sp h. We are able to substitute the first of them 
by considering 



d 2 H d 2 H 



dR 2 dLdR 



dL\ 



d 2 H 



(dxL) 



d XL dR\dR J sph 



(49) 



where we have omitted vanishing terms. Rearranging for 
(dL/di?) sp h and inserting into (48) we obtain 



dE\ 
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v3L/ \dR 2 

I d 2 H \ 
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m)] 


\dR 2 J 






J 



[(ail at)] o 



(50) 



The second part has been neglected as xl is supposed to 
be an almost conserved quantity even though radiation 
reaction is turned on. According to the analysis of [21] 
the corresponding error is of 3PN order. With [R]q in 
hand we can exploit the abovementioned quadratic ap- 
pearance of Pr within the Hamiltonian to obtain 



[Pr]o = 



[R]c 



1 dH 


2 


dH 




Pr 0Pr_ 


0,P R -+0 


[9(P R )_ 






(51) 



Hence, the algorithm gives [Pt]o and [Pr]o m a co- 
ordinate system where Lq (by definition) represents the 
z-axis. We omit the subscript "0" from here on as it is 
clear that we refer to initial data. Furthermore we use 
the labels R and D, denoting the binary separation, in- 
terchangeably 

The procedure just described is valid for both types 
of Hamiltonians. For the Taylor-expanded version we 
immediately obtain the final result in ADMTT gauge. 
Employing the EOB approach on the other hand is rather 
more involved. We shall briefly illustrate why. 

First, recall that in the EOB case the above algo- 
rithm is applied for the primed variables, which repre- 
sents a simplification as we discussed above. However, 
in the end we want to obtain the result in ADMTT co- 
ordinates. Moreover, we wish to specify the ADMTT 
distance, i.e. X, at which the initial momentum is to 
be determined. The canonical transformation ADM o 
EOB only mediates between the whole set of variables 
(X,P) -H> (X',P f ), thus providing the corresponding 
value of X' for use within the algorithm would already 
require the knowledge of the sought-after solution P. 
Hence, an additional step becomes necessary, for exam- 
ple a shooting method, to achieve a match with the given 
ADM distance. We take care of this for our actual initial 
data computations in order to compare numerical results 
at identical values of X. 

Translating EOB results to ADM coordinates involves 
a certain degree of inaccuracy because the transforma- 
tion is only known approximately (to 3PN order). Ta- 
ble|l]gives some conservative, rough error estimates based 
on PN convergence as shown in Fig. [T] The estimates 
already include the error contribution that would arise 
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Mass ratio 


Ro 1 M 




















6 


7 


8 


9 


10 


11 


12 


13 


14 


1:1 


3.1 


1.8 


1.1 


0.78 


0.55 


0.40 


0.30 


0.23 


0.18 


4:1 


2.2 


1.2 


0.74 


0.49 


0.34 


0.25 


0.19 


0.14 


0.11 


10:1 


0.99 


0.55 


0.34 


0.21 


0.15 


0.11 


0.08 


0.06 


0.05 



TABLE I: Estimate of the relative error in percent for EOB-based initial data Pt and Pr at various binary separations and 
mass ratios. 



from matching with a given ADM separation when ap- 
plying the shooting procedure. Apparently, the transfor- 
mation can be regarded more accurate for larger separa- 
tions and higher mass ratios. We conclude this section 
by presenting some examples for the orbital parameters 
Pt and Pr obtained by the algorithm. Figs. [2] and [3] 
show initial parameters for an equal-mass binary without 
spins. Significant deviations between the various combi- 
nations, that is combining one of the Hamiltonians, either 
Taylor-expanded (TH) or EOB Hamiltonian (EH), with 
one version of the flux function, either Taylor-expanded 
(TF) or Pade-resummed (PF), can be observed only at 
small separations. Some dependence on the choice of 
the flux function is observed in Fig. |3j although the re- 
sult depends more strongly on the choice of Hamiltonian. 
This effect becomes more pronounced for higher mass ra- 
tios, see Fig. [4j The plots clearly indicate the expected 
degradation of the TH parameters for small separations, 
which is caused by a z ero transition of the second deriva- 
tive d 2 H/dR 2 in Eq. (50), while the first derivative and 
thus d 2 H/ (dRdL) does not vanish. The resulting pole 
moves outwards to larger separations and spoils the data 
as the mass ratio increases. We would conclude that the 
accompanying drift in the Pt curves is also a result of 
the Taylor Hamiltonian's behavior. Clearly, it would not 
make sense to use TH data in the vicinity of the pole, 
Fig. |]d. 

The EOB data on the other hand appear sound and 



remain close to the P t values of [TO] (our Eq. (45)), but 
might be affected by the transformation error. In fact, 
in the equal-mass case there is concern that the trans- 
formation error exceeds the difference to the TH data. 
However, since the error goes down at higher mass ratios 
we expect the EOB data to be superior in these cases. 

We examine the data's behavior in the presence of 
spins at the end of Section [TV] 



IV. NUMERICAL RESULTS 

For a specific set of parameters, comprising masses, 
initial separation and spins, the initial data algorithm 
can give us tangential and radial momentum of a relative 
particle. The actual initial configuration, consisting of 
the individual locations and momenta of the black holes, 
can be obtained by placing the bodies appropriately in 
a coordinate frame by setting Pi = P-2 = P and e.g. 
letting the center of mass coincide with the origin. We 



are also free to apply rotations. 

As discussed the input parameters m a , X a , P a , S a are 
formulated in (almost) ADMTT gauge whereas the BAM 
code uses a different spatial metric 7y = ipQ 7^ with 
a conformally flat background metric 7^ = Sij on a 
maximal slice (K = 0) (so-called black hole puncture 
data |13j). In previous work it was found numerically 
that the two gauges are rather close to each other (e.g. 
[101 121)). However, it is by no means trivial nor en- 
tirely understood why PN-non-eccentric initial parame- 
ters should produce non-eccentric orbits in full GR. First, 
the process of translating the initial parameters to the 
actual initial data involves the solution of the Hamil- 
tonian constraint, which, for example, alters the mass 
parameters, cf. Table [V] Second, the modelled configu- 
ration lacks the gravitational-wave content which would 
be there if the binary had reached this state in a phys- 
ically realistic way. The constraint solution process and 
the ad-hoc specification of data at a given time result in 
artificial waves in the initial data, which e.g. affects the 
black-hole positions in the beginning of the simulation. 
Third, the gauge (lapse and shift) evolves in the starting 
phase, settling down to the final values after about 50M 
of coordinate time. These effects may not only outweigh 
the hard-won subtle differences between THTF, THPF, 
EHTF and EHPF initial data, but also make it difficult 
to reliably assess the quality of these data. 



Mass ratio 


ftmin,i/Mi 


r out /M 


1:1 


3/64 


774 


1:2 


3/56 


2067 


1:4 


1/96 


1236 



TABLE II: Grid setups of the different types of numerical 
simulations. The index i enumerates the black holes (i = 
1,2) and M = Mi + M2 is the total mass of the system. 
hminj is the resolution of the finest grid covering the i th black 
hole. r out is the position of the outer boundary. All results 
throughout the paper corresponding to the same mass ratio 
(but e.g. different initial separations or spins) were obtained 
with the setup given here. 



For the purpose of the latter we investigated several ap- 
proaches for measuring the eccentricity of the inspiralling 
binary's orbit. Following the ideas of O EU |43l [44], we 
employ two approaches based on the evolution of the sep- 
aration D(t) and the orbital frequency u>(t). The first 
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FIG. 2: (Color online) Tangential momentum obtained with the initial data algorithm for an equal-mass binary [M = 1, no 
spins). Also shown are reference functions "3PN circular" from [10] (cf. Eq. (45 1) and "PN evolved" from |21j . The right panel 
gives the deviation in percent with respect to the 3PN-circular curve. 
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FIG. 3: (Color online) Radial momentum obtained with the initial data algorithm for an equal-mass binary (M = 1, no spins). 
Compared are combinations of the two Hamiltonians (Taylor or EOB) with the two possible flux functions (Taylor or Pade). 
As expected, the THTF curve almost coincides with the reference function "PN evolved" from [21] . 



Mass ratio 


D/M 


Data type 


Pt / M 


P R 1 10" 3 M 


Interval /M 


Eccentricity 


1:1 


11 


THTF 


0.090110 


-0.71442 


300. 


.600 


0.002 






EHPF 


0.090266 


-0.71489 


300. 


.600 


0.003 


1:1 


7 


THTF 


0.123863 


-3.34124 


60. 


.120 


0.005 






THPF 




-3.13589 


60. 


. 120 


0.005 






EHTF 


0.125288 


-4.04062 


80. 


.150 


0.005 






EHPF 




-3.79255 


80. 


.150 


0.005 


2:1 


10 


THTF 


0.085615 


-0.79486 


250. 


.400 


0.002 






EHTF 


0.085753 


-0.80701 


250. 


.400 


0.003 


4:1 


10 


THTF 


0.061921 


-0.43393 


300. 


.700 


0.004 






EHPF 


0.061883 


-0.42106 


300. 


.700 


0.0025 



TABLE III: Simulations of black hole binaries without spin for equal and unequal masses. A local eccentricity measurement 
was performed in the time interval given. The uncertainty in the eccentricity is estimated to be ±0.001. 



consists of comparing the actual with an ideal (quasi- 
circular) inspiral curve D r {t) or w c (i), which is ob- 
tained by averaging over the eccentricity oscillations us- 



ing an appropriate fit function (also see Fig. [9j|. The 
resulting eccentricities are referred to as ejj and e u . 
Specifically, we employ a model function according to 
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D c {t) — J2t=i a i (*Af — *) 1 ^ 2 an< ^ eccentricity is com- 
puted as e D (t) = [D(t) - D c {t)]/D c {t) 21 . The time 
tM entering the formula for the ideal separation D c {t) 
is the merger time which in our case is estimated by 
the fitting procedure. The resulting value agrees very 
well with results from different methods for measuring 
the merger time, for example the time when the am- 
plitude of ^4 reaches its maximum. As a model func- 
tion for the frequency we use a fourth order polynomial, 
w c(i) = Si=o bit*) following [35]. In this case, the ec- 
centricity is calculated via e u (t) = [oj(t) — u) c (t)]/2u) c (t). 
The fit to the model function only works well if we ex- 
clude initial gauge adjustments and it fails at merger time 
(D c (t = thi) — in the D- method and the fourth or- 
der polynomial ceases to capture the behavior of the fre- 
quency in the w-method). We therefore have to choose 
an appropriate time interval for the averaging procedure. 
Figures [5] [8] show graphs for eccentricities varying with 
time resulting from the D-method according to the for- 



mula given above. The global extremum of each curve, 
usually located at early times, is then used to determine 
a time-independent eccentricity value of the data which 
will serve to compare different runs with each other. 

A second method tries to capture the oscillations them- 
selves by presuming some functional dependence, i.e. a si- 
nusoidal fluctuation on a certain background function like 
D(t) = Ax + A 2 t + B sin(£! i + <po), where A u A 2 , B, ft 
and 4>q are fit parameters. This works best for the time 
derivatives D{t) and tu(t) instead of the original func- 
tions. The resulting eccentricity can be estimated to 
be ejj(t) « — B cos(f2 < + (j)o)/QoD(t). The procedure 
is equally applicable to the derivative of the frequency 
curve, and leads to an eccentricity denoted by ecj{t). 
Since the assumed "capture" functions are only able 
to trace the curves D(t) and ui(t) locally (in a very 
short time interval), eccentricities computed with this 
approach turned out to be less useful for our purpose. 
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All in all, the different methods (e_o, e u , and e^) 
do not always give consistent results for very low eccen- 
tricities, for which the quality of the fit function becomes 
decisive. Eccentricity values given in the present paper 
are based on eu, with error estimates based on the vari- 
ations in the extrema of eo- 

Numerical evolutions were performed with the BAM 
code [ini ES] which uses the BSSN formulation of the 
Einstein equations [37J EE] an d the method of moving 
punctures (49J, [50] . Spatial derivatives are 6 th order accu- 
rate [5T] and a 4 th order Runge-Kutta scheme is employed 
for time stepping. We use "1+log" slicing for determin- 
ing the lapse function [52] and the "000" -version of the 
T-driver shift condition (53[ [54] . In order to resolve the 
black hole regions as well as the wave zone sufficiently 
well, mesh refinement with nested boxes is used. De- 
tails of the grid setups used for numerical simulations 
are given in Table [TTJ Note the choice of resolutions for 
the black holes: the number of refinement boxes is cho- 
sen for each black hole individually in order to obtain the 



A. Numerical results for different mass ratios and 
vanishing spin 

Table [HT] summarizes the outcome of eccentricity mea- 
surement for various simulations with equal and unequal 
masses without spin. The initial separation was partially 
chosen to allow direct comparison with (22 [25]. Fig. |5] 
shows the time dependence of for mass ratio 4:1, giv- 
ing a visual impression of the oscillations of the inspi- 
ralling orbit, while Fig. [6] shows results for 1:1. 

We draw the following conclusions from the numerical 
results. Our techniques for measuring eccentricity give 
compatible results in most cases, but in some cases the 
accuracy of ±0.001 is not quite sufficient to reliably re- 
solve the smallish differences. Unless the eccentricities 
deviate by amounts larger than the uncertainties we can- 
not name one of the data the better choice. In particular, 
we examined the consequences of the choice of flux for the 
smallest initial separation D — 7M case, where the dif- 
ferences should be largest. At the given level of accuracy 
we cannot discern a significant difference between the TF 
and PF methods. 

For initial separations of D = 10M and D = 11M, 
the eccentricity ranges from e = 0.002 to 0.004 for com- 
paratively small variations of about 1% in P t and 2% in 
Pr. This reflects the high sensitivity of the eccentricity 
to small changes in the momenta. The minimal eccen- 
tricity obtainable by the method is about e = 0.0025 for 
mass ratios from 1:1 to 4:1. 

For smaller initial binary separation the PN initial data 
become less appropriate since the underlying equations 
approach their limit of validity. Moreover, the deviation 
in the gauges, that is ADMTT gauge on the PN side ver- 
sus maximal slicing in the numerical simulation, grows. 
The numerical example for D = 7M shows that the ec- 
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FIG. 5: (Color online) Dynamics and eccentricity for mass 
ratio 4:1 (D — 10 M). Panel (a) shows the oscillating decrease 
of the binary separation for the THTF data. The darker part 
of the line indicates the time interval used for the fit required 
for e_o. In panel (b) the eccentricity is compared between 
Taylor-expanded and EOB-Hamiltonian data. 



centricity grows to e = 0.005 for equal masses, which is 
twice as large as for D = 11M. Considering how close 
D = 7M is to the merger, it is surprising how well the 
PN initial parameters work in this context. 

Going to higher mass ratios the putative advantage of 
EH data first seems to show at 4:1, although the dif- 
ference in initial momentum to TH data is quite small. 
Given the available data, the smaller eccentricity might 
be accidental, but we believe there is a systematic trend 
based on corresponding results for spinning binaries. 
(Numerical simulations at 10:1 have recently become 
available |55j . so the behavior for increasing mass ratios 
could now be investigated further.) 

Utilizing the same Hamiltonian, our TH data are ex- 
pected to give results similar to the PN-evolved data, 
and this is indeed the case. For example, [21] obtained 
e D = 0.002±0.001 for the 1:1 case starting at D = 11M. 
Our TH results are also consistent with the outcome of 
PN-evolved data at higher mass ratios, for example [25] 
found e = 0.003 for mass ratio 2:1 and e = 0.005 for the 
4:1 case at D = 10M. 
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FIG. 6: (Color online) Eccentricity for mass ratio 1:1 (D = 
11M). Compared are runs for PN-circular data, PN-evolved 
data, and the new THTF and EHPF data. The THTF run is 
very similar to the PN-evolved run as expected, while for 2PN 
circular data the eccentricity is up to a factor of five larger. 



In [21 , which introduced the PN evolution method to 
obtain initial parameters with reduced eccentricity, the 
PN evolutions were started with 3PN circular parame- 
ters of [10] . and the result was compared to numerical 
simulations starting from 2PN (one order lower) circu- 
lar parameters based on [22j E3] . The PN-evolved data 
showed an improvement by a factor of up to five in the 
eccentricity compared to the 2PN circular parameters for 
mass ratio 1:1. In Fig. [6] we compare the result for 1:1 
for 2PN and 3PN circular parameters, the PN-evolved 
parameters, and the new THTF and EHPF parameters. 
The THTF run is very similar to the PN-evolved run as 
already pointed out. As stated in Table III the EHPF 
parameters result in about 50% larger eccentricity (for 
equal masses). For numerical simulations starting from 
the 3PN instead of 2PN circular parameters the improve- 
ment is less, about a factor of two for mass ratio 4:1 
and 1:1 as shown in Figs. [5] and [6j Incidentally, note 
that the 2PN circular parameters are based on harmonic 
coordinates while the 3PN circular parameters refer to 
ADMTT. 



Our final comment on the vanishing spin case is that, 
as shown in Figs. [5] and |6j there still is a distinctive or- 
bital eccentricity for mass ratio 4:1, while for 1:1 the ec- 
centricity is reduced to a level where the oscillations are 
less clearly associated with the orbital motion. At least 
for 4:1 there seems to be room for improvement. 



B. Numerical results for equal masses and 
non-vanishing spin 

Let us now turn to spinning black hole binaries for 
equal masses. We would like to emphasize that as soon 
as spins are present eccentricity is inevitable. This is due 
to spin-spin interaction and can be seen from the PN 
equations of motion ^ , which predict a precession of the 
spins with a frequency determined by (18 1. This preces- 
sion influences the orbital motion via Eq. ( 16 ). Moreover, 



we find an even more direct impact on the orbital motion 
as a consequence of the i?-dependence in ( 18 1 which di- 
rectly affects the evolution of P in The latter mech- 
anism is also relevant when both spins are parallel to 
the orbital angular momentum. Despite the constancy 
of L = X x P in these cases (ignoring radiation for the 
time being) , the orbital variables X and P would still be 
allowed to exhibit the corresponding oscillations. 

Table HVI a illustrates this issue for an evolution with 
the conservative PN equations. (This is different from 
the PN-evolved parameters we discuss elsewhere which 
include flux terms.) The test case consists of an 
equal-mass binary with spin settings according to Si = 
Grrii/c 2 ■ \L. We applied suitable initial data from the 
above algorithm which would give circular orbits if only 
Hgo was involved (in particular we dropped the [Pr]o 
component). The observed radial fluctuations are of or- 
bital frequency (when taking periastron advance into ac- 
count) and can be explained by means of Eqns. ( 18 1- 
( 20 1 . Since neither L nor the spins precess in the cur- 



rent constellation, the spin-spin contribution just repre- 
sents a spherically symmetric deformation term H^s = 
const • R~ 3 . Due to its simplicity in this particular case 
the correction could in principle be taken into account 
when computing the initial parameters. For general spin 
configurations, however, this is not advisable, so we do 
not pursue this idea here but return to it later. Instead 
we can regard the numbers in Table |IV| a as an estimate 
of the minimal eccentricities to be achieved in fully nu- 
merical simulations. (For \ — 0: the eccentricity is ana- 
lytically zero, and the small deviations from zero in Ta- 
ble IV a are numerical error of the integration method.) 

The aligned and anti-aligned scenarios were explored 
with several runs. The main purpose was to find out 
which of the data type (THTF, THPF, EHTF, EHPF) 
leads to minimal eccentricity. The observation that the 
choice of the flux does not play an important role seems 
to carry over to the spinning case, so we focussed on the 
effect of the Hamiltonian. Starting evolutions at separa- 
tions between D = 10 . . . 12M we varied the spin magni- 
tude. Table [V] gives the corresponding initial parameters 
and also shows the effect of constraint solving. For spin- 
ning holes the bare masses fh\ and fri2 are reduced to keep 
the total energy constant. In our notation m a replaces 
the capital M a often used in the numerical literature in 
order to maintain the PN notation. 
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X 


en (PN evolution) 
D = WM D = 8M 


°) X 


Cd { i n) 




-0.25 


0.0020 


0.0040 


-0.25 


0.015 ± 0.002 


0.005 ±0.001 


0.00 


1 ■ 10" 7 


3 • 10" 7 


0.00 


0.002 ±0.001 


0.003 ±0.001 


0.25 


0.0015 


0.0024 


0.25 


0.015± 0.003 


0.0015 ±0.0005 


0.50 


0.0050 


0.0083 


0.50 


0.034 ±0.005 


0.005 ±0.001 


0.85 


0.014 


0.022 


0.85 


0.07 ±0.01 


0.014 ±0.003 



TABLE IV: a) Eccentricities from conservative PN evolution for aligned equal spins at D 
from corresponding simulations with BAM at comparable binary separations, cf. Fig. [7] 



10M and D = 8M. b) Results 



Mass ratio 


X 


D / M 


Data type 


Pt 1 M 


P R / 10~ 3 M 


mi 


mi 


mi 


m2 


1:1 


-0.25 


12 


THTF 


0.086803 


-0.60477 


0.5 


0.5 


0.4757 


0.4757 








EHPF 


0.086459 


-0.57913 


0.5 


0.5 


0.4757 


0.4757 


1:1 


0.00 


11 


THTF 


0.090110 


-0.71442 


0.5 


0.5 


0.4872 


0.4872 








EHPF 


0.090266 


-0.71489 


0.5 


0.5 


0.4872 


0.4872 


1:1 


0.25 


12 


THPF 


0.083322 


-0.47391 


0.5 


0.5 


0.4758 


0.4758 








EHPF 


0.083858 


-0.49460 


0.5 


0.5 


0.4758 


0.4758 


1:1 


0.50 


11 


THPF 


0.085921 


-0.54782 


0.5 


0.5 


0.4328 


0.4328 








EHPF 


0.087194 


-0.60170 


0.5 


0.5 


0.4328 


0.4328 


1:1 


0.85 


10 


THPF 


0.087364 


-0.60110 


0.5 


0.5 


0.2564 


0.2564 








EHPF 


0.090080 


-0.71349 


0.5 


0.5 


0.2563 


0.2563 


4:1 


0.00 


10 


THTF 


0.061921 


-0.43393 


2.0 


0.5 


1.9788 


0.4761 








EHPF 


0.061883 


-0.42106 


2.0 


0.5 


1.9788 


0.4761 



TABLE V: Initial parameters used in the simulations with spin for equal masses. The same quantities are shown for a 4:1 run 
for comparison. The total mass is M = mi ± m,2. 



Fig. J71 depicts the resulting eccentricity graphs while 
Table |IV| b provides the corresponding numbers. Fig. [8] 
shows additional details for the data in Fig. [7J including 
results obtained by other methods. 

The main observation is that the eccentricity increases 
significantly as x IS changed from zero towards positive 
and negative values. For TH data, the minimum is near 
X = 0, while for EH data the minimum is closer to 
X = 0.25 than x — 0- The EOB data yield remarkably 
low eccentricities, and for x ^ 0-25 the values are very 
close to the PN-limit estimate of Table EDa. The TH 
data perform poorly in comparison, with up to five times 
larger eccentricity for |vj > 0.25. It would be interesting 
to extend this study to more negative values of x- I n or ~ 
der to illustrate these results Fig. [9] shows the black- hole 
orbits of the simulation with x — 0.85. Strong deviations 
from a spiral can be observed for the Taylor-Hamiltonian 
data while the EOB data produce a rather well-formed 
inspiral. Note that the motion takes place in the z = 
plane. 

As a first step towards exploring more general spin 
configurations, we consider one particular example with 
arbitrarily chosen spin vectors in Fig. [7]f. In this case 
the motion does not remain in the z = plane and the 
orbital plane shows considerable precession, see Fig. [10] 
Both spins are of magnitude x = 0.5 with a smaller but 



positive component in the z-direction, and the eccentric- 
ity is indeed in the range that we would expect from the 
graphs Fig. [7]c-d where spins of magnitude 0.25 and 0.5 
are aligned to the z-direction. This indicates that our 
initial data are not restricted to aligned spins only, but 
also lead to small eccentricity for more general spins. 

Finally, let us look at the spin-spin interaction issue 
once more. As stated in the discussion of the aligned 
spin setup, the spin-spin interaction term simplifies in 
this case and might be incorporated in the initial data 
computation. This was done in Fig. [8] and indeed eccen- 
tricity is reduced by a factor two for the EH data. This 
approach, however, cannot amend the poor performance 
of the TH initial parameters. 

After noting the effect of the different initial data 
types, let us investigate how the large differences can be 
explained by having a closer look at the particular values 
of P t and Pr. In the equal-mass setup without spins the 
differences among the initial parameters were compara- 
tively small even up to close binary separations, Figs. [2] 
and [3] When incorporating spins (parallel to L) we ob- 
serve the modifications to tangential and radial momen- 
tum shown in Fig. 11 As a first finding the tangential 



component experiences significant changes all over the 
spin-parameter space. This can be explained by means 
of frame dragging. Spacetime is twisted in the vicinity 
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FIG. 7: (Color online) Performance of PN initial data with spin in numerical simulations. In the present equal-mass con- 
figurations (mi = m2 = 0.5), both spin vectors point in the same direction parallel to the orbital angular momentum and 
have the same length, i.e. Si = Gnii/c 2 ■ \L for panels (a) - (e). We compare the eccentricity e_o between TH and EH data 
for various values of X- In panel (f) we set up a more general spin configuration with Si = mfxi (0, — cos 45°, cos 45°) and 
S2 = rn\ X2 (— cos 60°, 0, cos 30°) with Xi ~ X2 = 0.5. In order to simplify comparisons of the plots, all axes are scaled the same 
way, and tM denotes the moment of merger. Large amplitudes in the beginning of the graphs are due to gauge adjustments in 
the code and should be ignored, while strong drifts near the end of some curves arise when the fit function (needed to compute 
e_o) becomes inappropriate. 
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FIG. 8: (Color online) Zoom-in plots for some of the cases of Fig. [7] comparing to eccentricities obtained with other methods, 
a) For equal masses and spin x — 0.25, 2PN-circular data (Pr = 0) amended with spin (see [56]) happen to give larger 
eccentricity than the EH data with properly chosen Pr. b) For equal masses and spin \ = 0.85, the benefits of including 
spin-spin interactions in the initial data computation are shown. 



of each spinning hole and the other object is forced to 
follow this twist. When both spins are aligned with the 
orbital angular momentum each black hole tries to ac- 
celerate the other on its orbit, with the effect that less 
orbital momentum P t is needed for a circular inspiral. 
(Recall that it is the velocity which is responsible for 
the orbit's shape and that there is a complicated, spin- 
dependent relation between X = dH(X, P, Si, S2)/dP 
and the momentum.) 

The TH and EH surfaces almost coincide in the non- 
spinning case, which marks the center of the plots. 
With increasing spins the Taylor-expanded data exhibit 
a stronger reaction than the EOB ones. Recall that com- 
paratively small changes in P t can have strong effects on 
the orbit. The TH and EH surfaces for Pr are almost 
identical for a certain part of the parameter space, in 
particular as long as both spins are positive. For strong 
negative spins, especially for the dominant spin of the 
bigger black hole at mass ratio 4:1, large changes in P^ 
occur. Obviously, the Taylor-expanded data are affected 
most, resulting in a behavior similar to Fig.|4]d, such that 
no reasonable solution is obtained beyond a certain spin 
threshold. Considering their poor performance in numer- 
ical runs, we conclude that TH data are as inappropriate 
in this regime as in the non-spinning case at high mass 
ratios. 

The EOB-based data do surprisingly well recalling the 
difficulties that arise from the transformation EOB — > 
ADM. The systematic simplifications we imposed, in par- 
ticular the approximation concerning the binary sepa- 
ration R' entering the spin-orbit Hamiltonian, Eq. ( 16 ) 



resp. (40), do not seem to significantly degrade their per- 



formance. 



V. CONCLUSIONS 

Based on the work of [25] we presented several varia- 
tions of an (almost) analytical algorithm capable of de- 
livering quasi-spherical initial data for black-hole binary 
systems. By employing two different PN Hamiltonians 
in combination with two versions of energy-flux functions 
we constructed four types of initial parameters. For the 
resulting initial data we performed numerical simulations 
within the moving puncture framework with the BAM 
code, both for non-spinning and spinning binaries, where 
we focussed on the class of aligned spins. An eccentric- 
ity measurement provided information on the quality of 
the data. In summary, the eccentricity of the initial data 
considered grows for smaller binary separations, larger 
mass ratios and larger spin magnitudes. 

In the non-spinning case, differences in the initial pa- 
rameters P t and Pr as well as in the resulting eccen- 
tricities were found to be comparatively small and often 
beyond the accuracy of our measurement. However, TH 
data happen to suffer from a parasitic pole at higher mass 
ratios. In particular due to the EH data we were able to 
meet or even to obtain slightly lower eccentricities than 
those resulting from the PN-evolution method of [21] . As 
in [21 , a non-zero radial component of the initial momen- 
tum is required for best results. 

As soon as spins are incorporated, PN theory shows 
that non-eccentric inspirals are not possible in general 
due to spin-spin interactions, but we demonstrated that 
EH initial data are suitable to minimize the radial oscil- 
lations to a certain degree, at least for special, aligned 
spin constellations. We considered one example for non- 
aligned spins, for which a similarly small eccentricity was 
obtained. Further runs with arbitrary spin directions are 
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-4 -2 2 4 500 1000 

X/M t / M 




-4 -2 2 4 500 1000 



X/M t/M 

FIG. 9: (Color online) Puncture tracks (left panels) and puncture distance (right panels) for the \ = 0.85 simulation from 
Fig. [7]e. The upper plots pertain to the Taylor-expanded data (THPF) and give a large eccentricity of eo = 0.07, whereas 
the EOB data (EHPF) yield much smaller fluctuations in the binary separation leading to eo = 0.014 (lower plots). The 
right-hand-side graphs also indicate the procedure of eccentricity measurement. The grey/lighter curve is the fit function D c (t) 
representing an ideal, non-eccentric inspiral. 




FIG. 10: (Color online) Orbital motion for the general spin example from Fig. [7]f. The orbital plane precesses and we observe 
a spin kick of the final black hole. 
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a) mass ratio 1:1 b) mass ratio 1 : 1 




c) mass ratio 4:1 d) mass ratio 4:1 




FIG. 11: (Color online) Modifications of PN initial parameters Pt and Pr for a binary at a fixed initial separation of 8M 
when spins are involved. The plots show THTF (blue/dark) and EHTF (red/light) initial data surfaces, for which the z- 
component of the spins is varied in the interval [— 1, +1], The zero-spin case marks the center of the plots, while e.g. the 
Xiz = X2 z = 1 case corresponds to maximally-spinning black holes with spins aligned to the orbital angular momentum. 
Changing to Pade-resummed flux does not cause visible modifications of the surfaces. 



planned for future research. As expected, the eccentric- 
ities in simulations with spins turned out to be larger 
than the spinless ones. The TH data gave comparatively 
poor results, which again can be considered an effect of 
the pole mentioned above. 

To end with a concrete suggestion for practical imple- 
mentations, a good starting point is the EOB Hamilto- 
nian with the Taylor expanded flux, especially for non- 
zero spin, but also for mass ratios 4:1 and larger. For 
smaller mass ratios the Taylor Hamiltonian gives smaller 
eccentricities, but the eccentricity for both EH and TH is 
quite small. The main complication in the EOB method 



is the need for the transformation to ADMTT coordi- 
nates for the numerical initial data. Table III gives 



explicit values for initial parameters for several non- 
spinning binary configurations that can be used directly 
in simulations. Use the combination of P t and Pr for 
which the eccentricity given in the last column of the 
table is lowest. Initial data for several binary configura- 
tions with (anti-)aligned spins can be read off table |v| 
For all cases with \ 7^ 0, the EHPF data leads to the 
lowest eccentricity. 

Further investigations suggest themselves. The PN- 
evolution method of 21 depends on the quality and char- 
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acteristics of the Hamiltonian evolution method. It is 
natural to revisit this topic with the insights gained in the 
present work for unequal masses and spins, and also to 
compare with the PN-evolution methods of [57] . Our 
study centers on the PN method of [21] and numerical 
evolutions. A refined method called "post-post-circular" 
initial data has been suggested in |57j . which is worth 
investigating in the context of numerical evolutions as 
well. 

Apart from evolutions, it could be interesting to com- 
pare to other diagnostics for the eccentricity of orbits, 

e.g. HUSOES]. 

An alternative to the TH and EH methods considered 
here is the semi-analytic puncture evolution approach to 
model inspiralling black holes [50] . It is natural to expect 
that the resulting orbital parameters are well adapted to 



numerical puncture evolutions, and we plan to carry out 
a quantitative study. 
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